#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <math.h>
#include <string.h>
#include <errno.h>
#include <time.h>
#include <sys/time.h>

#define NBR_THREADS 100
#define MY_RAND_MAX 2147483647

struct BoxMullerState
{
    double x1, x2, w, y1, y2;
    int useLast;
    struct drand48_data random;
};

void *func_thread(void *);
void *xmalloc(unsigned long int);
double max(double, double);
int old_gnu_rand_r(unsigned int *);
void initBoxMullerState(struct BoxMullerState*);
double boxMullerRandom(struct BoxMullerState*);

long int M;
double *trials;
double S, E, r, sigma, T;

int main()
{
    pthread_t threads[NBR_THREADS];

    int i;
    long int long_i;
    double mean = 0, std_dev = 0, conf_width, conf_min, conf_max;

    scanf("%lf %lf %lf %lf %lf %ld", &S, &E, &r, &sigma, &T, &M);

    trials = xmalloc(M * sizeof(double));

    for(i = 0; i < NBR_THREADS; ++i)
    {
        if(pthread_create(&threads[i], NULL, func_thread, (void *)i) < 0)
        {
            fprintf(stderr, "%s: %d: %s\n", __FILE__, __LINE__, strerror(errno));
            exit(EXIT_FAILURE);
        }
    }

    for(i = 0; i < NBR_THREADS; ++i)
    {
        pthread_join(threads[i], NULL);
    }

    for(long_i = 0; long_i < M; ++long_i) mean += trials[long_i];
    mean /= M;

    for(long_i = 0; long_i < M; ++long_i) std_dev += (trials[long_i] - mean) * (trials[long_i] - mean);
    std_dev = sqrt( (std_dev / (M - 1)) );

    conf_width = 1.96 * std_dev / sqrt(M);
    conf_min = mean - conf_width;
    conf_max = mean + conf_width;

    printf("S\t%lf\n", S);
    printf("E\t%lf\n", E);
    printf("r\t%lf\n", r);
    printf("sigma\t%lf\n", sigma);
    printf("T\t%lf\n", T);
    printf("M\t%ld\n", M);
    printf("Confidence interval: (%lf, %lf)\n", conf_min, conf_max);

    return 0;
}

void *func_thread(void *msg)
{
    long int long_i, size = M / NBR_THREADS;
    int i = (int) msg;
    int start = i * size, end = (i + 1) * size;

    //unsigned int seed = time(NULL);
    struct BoxMullerState state;
    initBoxMullerState(&state);

    for(long_i = start; long_i < end; ++long_i)
    {
        //double random_number = ( (double)old_gnu_rand_r(&seed) ) / MY_RAND_MAX;
        double random_number = fabs(fmod(boxMullerRandom(&state), (double)1));
        double t = S * exp( (r - 0.5 * sigma * sigma) * T + sigma * sqrt(T) * random_number );
        trials[long_i] = exp(-r * T) * max((t - E), 0);
    }

    return NULL;
}

void *xmalloc(unsigned long int size)
{
    void *ptr = NULL;
    ptr = malloc(size);
    if(!ptr)
    {
        fprintf(stderr, "%s: %d: %s\n", __FILE__, __LINE__, strerror(errno));
        exit(EXIT_FAILURE);
    }
    return ptr;
}

double max(double a, double b)
{
    return a >= b ? a : b;
}

/* Reentrant random function from POSIX.1c.
   Copyright (C) 1996, 1999, 2009 Free Software Foundation, Inc.
   This functin is part of the GNU C Library.
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. */
int old_gnu_rand_r(unsigned int *seed)
{
    unsigned int next = *seed;
    int result;

    next *= 1103515245;
    next += 12345;
    result = (unsigned int) (next / 65536) % 2048;

    next *= 1103515245;
    next += 12345;
    result <<= 10;
    result ^= (unsigned int) (next / 65536) % 1024;

    next *= 1103515245;
    next += 12345;
    result <<= 10;
    result ^= (unsigned int) (next / 65536) % 1024;

    *seed = next;

    return result;
}

void initBoxMullerState(struct BoxMullerState* state)
{
    state->random.__init = 0;
    state->useLast = 0;

    struct timeval now;
    gettimeofday(&now, NULL);
    state->random.__x[0] = now.tv_usec;
}

double boxMullerRandom(struct BoxMullerState* state)
{
    if (state->useLast)
    {
        state->y1 = state->y2;
        state->useLast = 0;
    }
    else
    {
        do
        {
            drand48_r(&state->random, &state->x1);
            state->x1 = 2.0 * state->x1 - 1.0;
            drand48_r(&state->random, &state->x2);
            state->x2 = 2.0 * state->x2 - 1.0;
            state->w = state->x1 * state->x1 + state->x2 * state->x2;
        }
        while (state->w >= 1.0);

        state->w = sqrt((-2.0 * log(state->w)) / state->w);
        state->y1 = state->x1 * state->w;
        state->y2 = state->x2 * state->w;
        state->useLast = 1;
    }

    return state->y1;
}
